Dynamics of stable viscous displacement in porous media 
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We investigate the stabilization mechanisms of the invasion front in two-dimensional drainage 
displacement in porous media by using a network simulator. We focus on the process when the 
front stabilizes due to the viscous forces in the liquids. We find that the capillary pressure difference 
between two different points along the front varies almost linearly as function of height separation in 
the direction of the displacement. The numerical results support arguments that differ from those 
suggested earlier for viscous stabilization. Our arguments are based upon the observation that 
nonwetting fluid flows in loopless strands (paths) and we conclude that earlier suggested theories 
are not suitable to drainage when nonwetting strands dominate the displacement process. We also 
show that the arguments might influence the scaling behavior between the front width and the 
injection rate and compare some of our results to experimental work. 



I. INTRODUCTION 

Immiscible displacement of one fluid by another fluid 
in porous media has important applications in a wide 
range of different technologies. Most often mentioned 
is hydrology and oil recovery. From a theoretical point 
of view, the displacement process is very complex and 
hard to describe in detail. Especially, much attention has 
been paid to the rich variety of displacement structures 
that is observed. The displacement structures are found 
to depend strongly on the fluid properties like viscosity, 
interfacial tension, fluid flow rate, and wettability [1-4]. 

In drainage the primary process is the displacement 
of a wetting fluid by a nonwetting fluid in porous media. 
Consider a two-dimensional (2D) horizontal displacement 
of a less viscous fluid by a more viscous fluid. At high 
injection rates the front developing between the invading 
and defending fluid, is known to stabilize [3]. In con- 
trast, at extremely low injection rate the invading fluid 
generates a growing cluster similar to the cluster formed 
by invasion percolation (IP) [5-8]. The displacement is 
now controlled solely by the capillary pressure, that is 
the pressure difference between the two fluids across a 
meniscus. 

In this paper we address the question of how the in- 
vasion front stabilizes when no gravity forces are present 
(2D horizontal displacement). To do this, we have de- 
veloped a network model that properly simulates the 
dynamics of the capillary pressures due to the menisci 
along the front as well as the viscous pressure buildup in 
the fluids. From the simulations we have calculated the 
capillary pressure difference AP^l between menisci along 
the front separated a distance A/i in the direction of the 
displacement. Also calculated, is the capillary pressure 
in the orthogonal direction APc±, that is the capillary 
pressure between menisci at same height above the in- 



let but separated a horizontal distance AZ (see Fig. 1). 
Simulations show that assuming a power law behavior 
AP^ii K Aft,'', our best estimate of the exponent for a 
wide range of injection rates and different fluid viscosi- 
ties is K = l.OiO.l. This is a surprising result because 
the viscous force fleld that stabilizes the front, is non ho- 
mogeneous due to trapping of wetting fluid behind the 
front and to the fractal behavior of the front structure. 

We also presents arguments being supported by the nu- 
merical evidence that k ~ 1.0. The arguments are based 
upon the observation that nonwetting fluid displaces wet- 
ting fluid through loopless strands (see Fig. 9). As a con- 
sequence, we flnd that existing theories [9-12] not consid- 
ering this effect, are not compatible with drainage when 
nonwetting strands dominate the displacement process. 
We also conjecture that the result k ~ 1.0 may influence 
the scaling between the saturated front width Wg and the 
capillary number Ca- The capillary number is the ratio 
between viscous and capillary forces and in the following 
Ca = QpLnw/'^l- Here Q is the injection rate, E is the 
cross section of the inlet, and ^„u, is the viscosity of the 
nonwetting fluid. 

The effect of gravity on the front when the flu- 
ids have different densities has been thoroughly dis- 
cussed [13,9,14,15] and in slow drainage it is found that 
gravity may stabilize the front. Gravity causes a hy- 
drostatic pressure gradient in the ffuids and consider- 
ing a heavy nonwetting fluid below displacing vertically 
upwards a light wetting fluid, this gradient will stabi- 
lize the front. The displacement process corresponds ex- 
actly to IP with a stabilizing gradient [9,14,16] and the 
saturated front width Wg, has been shown to scale like 
Ws oc Bo~^ ^"^ . Here u is the correlation length expo- 
nent in percolation and Bo is the bond number indicating 
the ratio between gravity and capillary forces. 

A similar consensus concerning the stabilization mech- 
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FIG. 1. A schematic picture of the front that travels across 

the system from the inlet to the outlet. In the figure, APc|i 

is the capillary pressure difference between a meniscus at A 

and a meniscus at B separated a vertical distance A/i. In the 

orthogonal direction we calculate APc±, that is the capillary 

pressure difference between a meniscus at B and a meniscus 

at C separated a horizontal distance AL AP„u, and APi^, 

denote the viscous pressure drop going from A to B in the 

nonwetting and wetting phase, respectively. 



Aft, in the direction of the displacement. Here, d^ is the 
EucUdean dimension of the space in which the front is 
embedded (in our case d^ = 2) and Ah is assumed to 
be sufficiently large for scaling to be acceptable and less 
than Ws- They also argued that the pressure drop in the 
wetting phase AP^, must be linearly dependent on Ah, 
since the displaced phase is compact. In [11] Blunt et al. 
also suggested a scaling relation for AP„^, however, in 
3D they found AP^w oc Aft*/''+i. This deviates from the 
result of Xu et al. when ^e = 3. 

The paper is organized as follows. In Sec. II we de- 
scribe the network model used in the simulations. Sec. Ill 
contains the simulation results of APc|| and APc±, sup- 
porting the arguments we present in Sec. IV. In Sec. V 
we compare our findings to some experimental data and 
the conclusions are drawn in Sec. VI. At the end we 
have put Appendix A where we deduce the scaling re- 
lation between Ws and Ca using the ideas in [9] when 
not considering the effect of nonwetting fluid flowing in 
strands. 



anisms when viscous forces replace gravity forces has not 
yet been reached. In the literature the displacement has 
been related to IP [9,11,12], however, the scenario is more 
complicated than in the gravity case. Gravity is a uni- 
form force acting on the whole system, while the viscous 
force is local and fluctuates due to permeability varia- 
tions and fluid trapping in the porous medium. One 
standard approach is to separate the displacement struc- 
ture into two parts. One consisting of the frontal region, 
and the other consisting of the static structure behind. 
The frontal region of extent Ws , is assumed to behave as 
the spanning cluster in percolation. Consequently, it is 

assigned the permeability k x Wg , where t is the con- 
ductivity exponent in percolation. By applying Darcy's 
law and assuming that the stabilized front reaches a 
traveling-wave state according to Buckley-Leverett dis- 
placement [17], the scaling of the front width is found to 
behave as Ws oc Cq"". In the literature there exists two 
slightly different expression for a. In 3D Wilkinson [9] 
found a = i'/{l + t — (3 + v) where trapping of wetting 
fluid is assumed to be less important. Here (3 is the or- 
der parameter exponent in percolation. Later, Blunt et 
al. [11] suggested in 3D that a = i^/(l + 1 + u) which is 
identical to the result of Lenormand [10] discussing limits 
of fractal patterns between capillary flngering and stable 
displacement in 2D porous media. In Appendix A we 
present a simple method giving a = h'/{l + t — P + v) by 
applying percolation concepts on the frontal region when 
not considering that nonwetting fluid flows in strands. 

Recently, Xu et al. [12] used Wilkinson's arguments 
and deduced a scaling relation for the viscous pressure 
drops in the frontal region. They proposed that the non- 
wetting pressure drop AP„^ in the front (see Fig. 1) 
should scale as AP„,„ ex A/i*/'^+'''=^^^^'^/'' over a distance 



II. NETWORK MODEL 

The network model has been presented else- 
where [18,19] and therefore only its main features will 
be given here. 

In the simulations we have constructed the porous 
medium in two different ways. In the flrst way the porous 
medium is represented by a square lattice of tubes ori- 
ented at 45°. The tubes are cylindrical with length d. 
Each tube between the ith and the jth node in the lattice 
is assigned an average radius rij which is chosen at ran- 
dom in the interval [Aid, A2rf], where < Ai < A2 < 1. 
The randomness of the radii represents the disorder in 
the system. In the following this system will be referred 
to as the random radii lattice. 

In the second way the porous medium is constructed 
upon a square lattice inclined 45° where the distance be- 
tween each intersection in the lattice is of unit length. 
Around each intersection we draw a circle of radius A. 
To avoid overlapping circles the given A must be in the 
interval < A < 1/2. A node is placed at random in- 
side each of the circles and the nodes inside the nearest 
neighbor circles are connected by cylindrical tubes. Thus, 
as for the random radii lattice, four tubes meet at each 
node. We let dij denote the length of the tube between 
the ith and jth node, and the corresponding radius r^ 
is defined as r^ = dij /2a. Here a is the aspect ratio be- 
tween the tube length and the radius. In the simulations 
a = 1.25, hence, the tubes are 25% longer than they are 
wide. In this lattice the position of the nodes represent 
the disorder in the system, and therefore we will refer to 
it as the random node lattice. 

While every pair of nearest neighbor nodes are sep- 
arated an equal distance in the random radii lattice, 
the distance between two nearest neighbor nodes vary in 



the random node lattice. Especially, the shortest length 
scale, that is the minimum distance between two neigh- 
boring nodes, is less in the random node lattice. Conse- 
quently, we are able to generate more narrow fronts at 
higher injection rates in the random node lattice, than 
what is possible with the random radii lattice. Therefore 
the random node lattice is preferred at high injection 
rates where a flat front is generated. 

In both lattices the tubes represent the volume of both 
pores and throats, and there is no volume assigned to the 
nodes. The liquids flow from the bottom to the top of 
the lattice, and we implement periodic boundary condi- 
tions in the horizontal direction. The pressure difference 
between the bottom row and the top row defines the pres- 
sure across the lattice. Initially, the system is filled with 
a wetting fluid with viscosity /i^. The injected fluid is 
nonwetting and has viscosity fXnw > Mm- The viscosity 
ratio M, is defined as M = ^nw/ ^^w■ 

The capillary pressure Pc between the nonwetting and 
wetting fluid in a tube is given by Young-Laplace law 



Pc 



1 



I 



Ri R2 



(1) 



where Ri and i?2 are the principal radii of curvature of 
the interface (a meniscus) and 7 is the interfacial tension. 
In a cylindrical tube of radius r where i?i — R2, Eq. (1) 
reduces to Pc = (27/r) cos 9. Here 9 denotes the wetting 
angle between the nonwetting and wetting phases, and 
in drainage 9 is in the interval (0,7r/2). 

In the network model we treat the tubes as if they were 
hourglass shaped with effective radii following a smooth 
function. Hence, we let the capillary pressure become a 
function of the meniscus position in the tube and assume 
the Young-Laplace law (1) takes the form 



Pc^yil- cos(2^f )] 



(2) 



Here < x < d is the position of the meniscus in the tube 
where d is the tube length. We assume perfect wetting, 
i.e. 61 = 0. 

By letting pc vary according to (2), we include the 
effect of burst dynamics into the model [18]. This is 
particularly seen at low injection rates where the inva- 
sion of nonwetting fluid takes place in bursts accompa- 
nied by sudden negative jumps in the pressure (Haines 
jumps) [20-22]. The detailed modelling of the capillary 
pressure costs computation time. However, it is neces- 
sary in order to properly simulate the pressure behavior 
along the front. 

The volume flux qij through a tube from the ith to 
the jth node is found from the Washburn equation for 
capillary flow [23] 



Qij 



^ij '^ij -L 



fJ-ij 



(Apij - PcA 



(3) 



Here fc^ is the permeability of the tube (r? /8) and ct^ 
is the cross section (7rr? ) of the tube, /i^ is the effective 



viscosity given by the sum of the volume fractions of each 
fluid inside the tube multiplied by their respective viscosi- 
ties. The pressure drop across the tube is Apij = pj —pi, 
where pi and pj is the nodal pressures at node i and j 
respectively. The capillary pressure Pc,y is the sum of 
the capillary pressures of the menisci (given by Eq. (2)) 
inside the tube. A tube partially filled with both liquids, 
is allowed to contain either one or two menisci. For a 
tube without menisci Pc^ij = 0, and Eq. (3) reduces to 
that describing Hagen-Poiseuille flow with /ijj = /ii or 

We assume conservation of volume flux at each node 
giving 



Yl *J = °- 



(4) 



The summation on j runs over the nearest neighbor nodes 
to the zth node while i runs over all nodes that do not 
belong to the top or bottom rows, that is, the internal 
nodes. 

Eqs. (3) and (4) constitute a set of linear equations 
which arc to be solved for the nodal pressures pi, with 
the constraint that the pressures at the nodes belonging 
to the upper and lower rows are kept fixed. The set 
of equations is solved by using the Conjugate Gradient 
method [24]. 

During every simulation we held the injection rate Q 
fixed and calculate a time dependent pressure AP across 
the system. See Refs. [18,19] for details on how AP and 
the corresponding p^'s are found. 

Having found the pi's we calculate the volume fluxes, 
Qij, through every tube in the network, using Eq. (3). 
According to the g^'s we define a time step At, such 
every meniscus is allowed to travel at most a maximum 
step length Ax^axi during that time step. The menisci 
are then moved a distance {qij /aij)At and the pressure 
AP and the time lapse are recorded, before the pi's are 
solved for the new fluid configuration. Menisci that are 
moved out of a tube during a time step are spread into 
neighbor tubes. For details about how the menisci is 
moved into neighbor tubes see Refs. [18,19]. 

Numerical simulations show that Axmax must be of or- 
der 0.1 to calculate the variation in the capillary pressure 
when a meniscus travel through a tube. In all simulations 
presented here Axmax — 0.1, resulting in at least ten 
time steps to invade one tube with nonwetting fluid. This 
causes the computation time to increase dramatically and 
one displacement simulation on lattices of sizes presented 
in this paper takes typically between 3-15 hours on a 400 
MHz Pentium II processor. 



III. SIMULATIONS 

We have run drainage simulations at different injection 
rates and fluid viscosities to study the capillary pressure 



variations along the invasion front. Due to the huge com- 
putational effort that is necessary, the simulations have 
been limited to lattices of size 25 x 35 and 40 x 60 nodes 
(Sec. Ill A). We have also run some simulations where 
the lattice initially was filled with nonwetting and wet- 
ting fluid according to patterns which were generated by 
an IP algorithm (Sec. IIIC). In this way, we were able 
to study the capillary pressure along invasion fronts on 
lattices of 200 x 300 nodes. 



In 



every 



simulation, APc||, ^Pi 



and the front 



width between the invading and the defending fluid, was 
recorded. The front was detected by running a Hoshen- 
Kopelman algorithm [25] on the lattice and recognized 
as the set of tubes that contain a front meniscus be- 
tween the nonwetting and wetting phase, that is the front 
tubes. The front width w is defined as the standard 
deviation of the vertical distances between front tubes 
and the average position of the front. Let hi denote 
the vertical distances of the front tubes above the in- 
let, where i — \,...,nf and Uf is the total number of 
front tubes. Then at a particular time, we calculate 
w = [(1/n/) ^j(/ii — h)"^]^/"^, where h is the average of 
the hiS. 

AP^ii and APc± is calculated as follows. Consider two 
front menisci denoted by m and n with height hm and 
hi above the inlet (bottom row) at a distance Im and Z„ 
from the left boundary of the lattice. Assume that hm > 
hn, then we calculate the difference APc"^"{Ah,Al) — 
Pc ~ pT where Ah — hm — hn and Al = \lm ~ ln\- If 
instead /i„ > h^, we compute APc"™(Aft,, A^) = p^-p" 
where Ah = hn — hm- We only consider the front tubes 
containing one meniscus and where the nonwetting fluid 
invades the tube from below. Note also, that we always 
take the capillary pressure of the meniscus closest to the 
inlet minus the capillary pressure of the meniscus closest 
to the outlet. From above, we define AP^i as function of 
Ah as the average of APc™" over all pairs mn separated a 
distance Ah but different Al, i.e APc|| = (APc™"(Aft, = 
const. , At)). 

The capillary pressure difference in the orthogonal di- 
rection, APc_L, (parallel to the inlet) as function of Al 
is defined as the average of jAPg™"! over all pairs mn 
with equal height (Ah = 0) above the inlet when Al 
is held constant. Thus, in the above notation APt,_L = 
(|APe'""(0,AZ = const.)!). 

The simulations were performed with parameters as 
close as possible to experiments performed in [26]. In 
the random radii lattice we set the length d, of all tubes 
equal to 1 mm and the radii r of the tubes were randomly 
chosen in the interval O.OBd < r < d. In the lattices 
with random nodes we chose the positions of the nodes 
such that the length of the tubes were inside the interval 
0.2 < d < 1.8 mm. This gave us the radii of the tubes, 
defined by r = d/2a, where a = 1.25. For both types of 
lattices the interfacial tension was set to 7 = 30 dyn/cm, 
and the fluid viscosities were 0.10 P, 0.50 P or 10 P. 



TABLE I. Simulations performed on the random radii lat- 
tice of size 25 x 35 nodes and M — 100 {fi„m = 10 P, /i^j = 0. 10 
P). The table contains the number of runs at each Q and Ca 
and the calculated Ws- 



Runs 


Q 

(cm^/min) 


Ca 


Ws 


30 


0.050 


3.7x10-* 


5.5 ±0.5 


30 


0.10 


7.3 xlO-* 


4.3 ±0.4 


30 


0.20 


1.5x10"^ 


3.7 ±0.4 


30 


0.50 


3.7x10^^ 


3.0 ±0.3 


30 


0.80 


5.8x10"^ 


2.5 ±0.3 


30 


1.5 


1.1x10"^ 


2.4 ±0.2 



TABLE IL Simulations performed on the random node lat- 
tice of size 25x35 nodes and M = 100 {^„m = 10 P, /i^, = 0.10 
P). The table contains the number of runs at each Q and Ca, 
and the calculated Ws- 



Runs 



Q 

(cm^/n 



Ca 



10 
20 
20 
20 
20 
15 
15 
10 
10 



0.010 

0.030 

0.050 

0.10 

0.30 

0.50 

1.0 

2.0 

4.0 



1.0x10"** 
3.1x10"** 
5.2x10"** 
1.0x10"^ 
3.1x10"^ 
5.2x10"^ 
1.0x10"^ 
2.1x10"^ 
4.2x10^2 



4.3 ±0.6 
2.9 ±0.3 
2.5 ±0.2 

2.1 ±0.2 

1.4 ±0.1 

1.2 ±0.1 
0.9 ±0.1 
0.8 ±0.1 
0.8 ±0.1 



A. Capillary pressure behavior 

We have performed two series of simulations with vis- 
cosity ratio M = 100 and one series of viscosity matched 
fluids, M ~ 1. In all series the capillary number Ca, was 
systematically varied by changing the injection rate Q. 
Tables I, II, and III list Q, Ca and the type of lattice 
(random radii or random nodes) used in the different se- 
ries. Also shown are the calculated front width Ws, and 
the number of different runs we did at each Q to obtain 
reliable average quantities. 



TABLE III. Simulations performed on the random node 
lattice of size 40 x 60 nodes and M — 1 (unw ~ l-tw ~ 0.50 
P). The table contains the number of runs at each Q and Ca 
and the calculated w.. 



Runs 


Q 
(cm'^/min) 


Ca 


Ws 


10 


0.050 


1.6x10"'' 


7.5 ±1.5 


10 


0.10 


3.2x10"^ 


6.9 ±1.2 


15 


0.30 


9.7x10"^ 


5.2 ±0.5 


15 


0.60 


1.9x10"** 


4.4 ±0.5 


20 


1.2 


3.9x10"** 


3.8 ±0.5 


20 


2.4 


7.8x10"** 


3.0 ±0.2 


20 


4.8 


1.6x10"^ 


2.4 ±0.2 
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FIG. 2. AP^ii as function of Ah for some Ca with M = 100 
(Table I) and M = 1 (Table III). AP^y is the average of 
the different runs performed at each Ca, and the error bars 
denote the standard error of the mean. Inset: logiQ{APc\\) as 
function of logi(,(A/i) for C„ = 1.1 x 10"^ and Ca = 3.7 x 10^'' 
with M — 100. The solid lines were fitted to the curves and 
their slopes are given by k. 



Fig. 2 shows the calculated capillary pressure difference 
APc|| , in the direction of the displacement as function of 
height separation Ah. We have plotted the result for 
some of the simulations performed on the random radii 
lattice of 25 x 35 nodes with M = 100 (filled symbols) and 
for some of the random node lattice of 40 x 60 nodes with 
M = 1 (open symbols). In the inset of Fig. 2 the results 
for highest and lowest Ca with M — 100, are plotted in 
a logarithmic plot and fitted to straight lines. Assuming 
a power law behavior, wc find that at Ca = 3.7x10^"' 
and M = 100, AP^n oc Aft" and k = 1.0. The expo- 
nent K seems to decrease systematically with increasing 
injection rate, and at Ca ~ 1-1 x 10~^ and M — 100 our 
best estimate is k = 0.8. Similar results was found from 
the simulations performed with viscosity matched fluids 
(M = 1). The data points corresponding to Ah < 1 
tube length, is omitted in the calculations of the expo- 
nent in Fig. 2. At short distances we expect uncertainties 
in the result because of the finite length of the tubes in 
the lattice. 

In Fig. 2 we observe that AP^l increases more rapidly 
as function of Ah at high injection rates compared to 
the results at low injection rates. In the plot the effect is 
most significant when M = 100. At extremely low injec- 
tion rate we expect AP^l in Fig. 2 to approach zero and 
become independent of Ah. In this limit the capillary 
pressure of the menisci along the front are almost equal 
(capillary equilibrium). As seen from Fig. 2, we have not 
performed simulations with that low injection rate. In- 
stead the lowest Ca for M ~ 100 and M — I, corresponds 
to the injection rate where no clear stabilization of the 
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FIG. 3. The average position {x), of the front menisci in- 
side the tubes as function of the menisci's height h relative 
to the of bottom height of the front hmin. The plot shows 
result from simulations in Table II at Ca = 1.0x10" (solid 
line), 1.0x10"^ (dotted line), and 1.0x10"^ (dashed line). 
Inset: The corresponding APc± as function of AL The lat- 
tice size was 25 x 35, giving a maximum horizontal distance 
Al — 12.5 due to the periodic boundary conditions in the hor- 
izontal direction. The error bars denote the standard error in 
the mean. 



front was found due to the finite size of the system. 

At higher injection rates the viscous gradient stabilizes 
the front. The gradient results the capillary pressure of 
the menisci closest to the inlet to exceed the capillary 
pressure of the menisci further down the stream. This 
is indicated in Fig. 3, showing the average position (x) 
of the front menisci inside the tubes as function of their 
vertical height h, relative to the bottom height of the 
front ft-niin- (x) is plotted for high, intermediate, and low 
Ca for the simulations listed in Table II. From the figure 
we observe that at high Ca = 1.0x10"^ (dashed line), 
the menisci near ft-min is placed closer to the middle of 
the tube compared to the menisci ahead. Consequently, 
the capillary pressure of the menisci near ft-min will more 
likely be larger than the capillary pressure of the menisci 
away from hmin and therefore tubes near hmin will more 
easily be invaded. This will eventually stabilize the front. 
Remember that the tubes are hourglass shaped and most 
narrow at a; = 0.5 (see Eq. (2)). At low injection rate, 
Ca — 1.0x10^^ (solid fine), we approach the regime of 
capillary equilibrium giving almost no difference in (x) 
as function of ft — ftmin- 

For the three Ca's in Fig. 3 we have also calculated 
the capillary pressure difference in the orthogonal direc- 
tion, APc±, as function of horizontal distance, AL The 
result is shown in the inset of Fig. 3. Here we interpret 
APc_L as the horizontal correlations in the capillary pres- 
sure between menisci at same height. Recall that APcj^ 
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FIG. 4. logj^g(APc||(™s)) as function of logjQ(Ca) for the 
simulations performed on the random node lattice with 
M = 100 (top) and M = 1 (bottom). The slope of the solid 
line in the upper figure is 0.15. The error bars denote the 
standard error in the mean. 



contains terms like |p™ — p" | = \J{p''c ~P?)^) where p™ 
and p" denote the capillary pressure of two front menisci 
m and n, respectively. From the inset of Fig. 3 we see 
that at low Ca = 1.0x10"'* (solid line) the capillary 
pressure of two menisci at same height and a distance 
A/ < 7 apart, are correlated to each other because AP^i 
as not yet reached the constant plateau (A/ > 7) where 
the capillary pressures becomes uncorrelated. At short 
distances APc± approaches zero, indicating that neigh- 
boring menisci have equal capillary pressures. At high 
Ca = 1-0 X 10~^ (dashed line), we observe that the corre- 
lations are very short. Already for A/ > 1, APc± reaches 
the plateau and the capillary pressures of the menisci do 
no longer interfere. Thus, if we consider a narrow and a 
wide tube at same height, the viscous forces are strong 
enough to push the nonwetting fluid through both the 
narrow and the wide tube simultaneously. As a result 
nonwetting fluid will invade simultaneously everywhere 
along the front. Similar behavior is observed in the other 
simulations listed in Tables I and III at high Ca 



B. Effect of viscosity ratio on the capillary pressure 

Fig. 4 shows a log-log plot of APc|| taken at A/i = Ws, 
as function of Ca for the simulations performed on the 
random node lattice with M — 100 (Table II) and 



M = 1 (Table III). In the following AP^n at w^ is de- 
noted as AP^ii {ws)- If we ignore the effect of nonwetting 
strands and use the result presented in Appendix A on 
our problem, we have that APc||(ws) oc CaWg'^ by set- 
ting Ah = Ws in Eq. (Al). Here Wg oc Ca""" where 
a — v/{l + t — (3 + v) and n = t/v + 1 — (3/v ac- 
cording to Appendix A. By combining the two power 
laws we obtain APc||(ws) oc c^i/(i+t-/3+'') giving in 2D, 
AP,||(ii;,)ocC„°-29. 

If we assume a power law behavior between AP^n {wg) 
and Cq, our best result for the exponent is 0.15 ± 0.05 
when M — 100 in Fig. 4. Note that there seems to be an 
upper cut off at Ca ^ 1.0x10"^ where APc||(ws) stops 
growing. At Ca ^ 1.0x10"^ the front is typically flat 
and we approach the minimum width due to the finite 
length of the tubes (see Table II). In this limit we expect 
a cross over to another type of behavior. 

If it is difficult to confirm any power law when M — 
100, the result of M = 1 in Fig. 4 does not show any 
scaling behavior. Already for Ca ^ 1x10^^, APc|[(ws) 
reaches a plateau or even decreases. To explain the differ- 
ent behavior of A.Pc\\{ws) when M = 1 and 100, we first 
look at the strength of the capillary pressure drop across 
the front and second we compare that to the magnitude 
of APcj^ as function of Ca- 

To study the capillary pressure drop we have calcu- 
lated the average capillary pressure {Pc) in the frontal 
region as function of the relative height from the bottom 
of the front, [h — /imin)/ws. The height is normalized by 
dividing with the saturated front width Wg ■ In the simu- 
lations {Pc) was computed by taking the average of the 
capillary pressures of the front menisci at same height, 
h, above the inlet. Fig. 5 shows the result for two simu- 
lations with almost equal Ca but different M . One with 
M = 1 and Ca = 1.6 xlO"^ (Table III) and the other 
with M = 100 and Ca = 1.0 x 10"^ (Table II). If we con- 
sider the middle part of the front between the two vertical 
dashed lines in Fig. 5, we observe that the capillary pres- 
sure drop, ~Wsd{Pc)/dh, over a length Ws in the front, is 
higher for M — 100 than for M — 1, even though the cap- 
illary numbers are almost equal. In both simulations a 
typical narrow front with a compact displacement struc- 
ture developed. On average, —Wsd{Pc)/dh must equal 
the difference between the pressure drops taken in the 
nonwetting and wetting part of the front over a length 
Ws (see Fig. 1). When the nonwetting and wetting fluid 
have equal viscosities the pressure drops in the nonwet- 
ting and wetting part of the front is about the same, ex- 
plaining the smaller capillary pressure drop when M = 1 
than when M = 100 in Fig. 5. 

Let us now study the behavior of APc^. Simulations 
show that APcj^ as function of Al docs not change much 
when comparing simulations performed at equal Ca with 
M — 1 and M = 100. Especially, the constant plateau 
where the capillary pressures are uncorrelated (see inset 
of Fig. 3), has the same value. This is illustrated in 
Fig. 6 where we have plotted the plateau of APcj^ versus 



2400 



3.2 



C,=1.6xlO ,M=1 

--- C =1.0x10"', M=100 




(h-h . )/w 

FIG. 5. (Pc) in the frontal region as function of the relative 
height from the bottom of the front. The height distance is 
normalized by dividing with the saturated front width Ws- 
The vertical dashed lines indicate the region where (Pc) is 
approximately linear. The error bars denote the standard 
error of the mean. 
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FIG. 6. The logarithm of the plateau of APc± versus the 
logarithm of Ca for M — 100 (circles) and M — 1 (boxes) 
corresponding to simulations listed in Table II and Table III, 
respectively. The slope of the solid line is 0.2. See also the 
inset of Fig. 3. 



Ca in a logarithmic plot for simulations with M — 100 
(Table II) and M = 1 (Table III). From the figure we 
observe that the plateau does not depend on M. As a 
side mark, we notice that there seems to be a power law 
between the plateau of APc± and Ca, which we indicate 
by the straight line in Fig. 6. The slope of the line is 0.2. 
From the above discussion we draw the following con- 
clusion. Consider two parallel and horizontal lines inter- 
secting the front and let the lines be separated a verti- 
cal distance Wg- When M = 1 we have found that the 
capillary pressure drop between the lines is small due to 
the equal fluid viscosities (Fig. 5). However, the magni- 
tude (plateau) of APci_ , is found to be the same as when 
M = 100 (Fig. 6). Thus, when M = 1 the relative small 
capillary pressure drop is annihilated by the magnitude of 
the capillary variations in the horizontal direction, APc±_ . 
This destroys a possible power law behavior of AP^i (ws) 
when M = 1 in Fig. 4. When M = 100, the capillary 
variations are too small to annihilate the larger capil- 
lary pressure drop there, giving the increasing function 
APj,||(ws). If we divide the capillary pressure drop, cal- 
culated in Fig. 5, with the plateau of APc± in Fig. 6, we 
find that the ratio is a factor three lower for M = 1 than 
for M = 100 at Ca ^ 1.0 x lO"-"^. 



C. Capillary pressure on IP patterns 

We have studied the capillary pressure along the front 
of patterns generated by an IP algorithm with a stabiliz- 
ing gradient. The patterns were loaded into our network 
model, and the simulations were started from that point. 
Using this method, we were able to perform displacement 
simulations in a short period of time on patterns gener- 
ated on lattices of 200 x 300 nodes. The result of these 



simulations are based on the assumption that the gen- 
erated patterns are statistically equal to the structures 
that would have been obtained in a corresponding com- 
plete displacement simulation. 

The IP algorithm was performed on the bonds in a 
square lattice with the bonds oriented at 45°. Hence, the 
bonds correspond to the tubes in our network model and 
an occupied bond refers to a tube filled with nonwetting 
fiuid. Each bond were assign a random number /^ in the 
interval [0, 1] where ij denote the bond between the ith 
and the jth node in the lattice. A stabilizing gradient g 
was applied on the lattice giving an occupation threshold 
tij of every bond like, Uj — fij + ghij [9,14]. Here hij 
denotes the height of bond ij above the bottom row. The 
occupation of bonds started at the bottom row, and new 
bonds were occupied until the invasion front reached the 
top row. There was periodic boundary conditions in the 
horizontal direction. The next bond to be occupied was 
defined as the bond with the lowest threshold value from 
the set of empty bonds along the invasion front. The 
invasion front was found by running a Hoshen-Kopelman 
algorithm on the lattice. 

We generated four IP patterns with g = 0.05 and dif- 
ferent sets of random numbers fij. When the invasion 
front became well developed with trapped (wetting) clus- 
ters of all sizes between the size of the bonds and the 
front width, the structures were loaded into our network 
model. Fig. 7 shows one of the generated IP patterns. 

The loading was performed by filling the tubes in the 
network model with nonwetting and wetting fluid ac- 
cording to occupied and empty bonds in the IP lattice. 
Furthermore, the radii r^ of the tubes were mapped 
to the random numbers fij of the bonds like, r^ = 
[Ai-|-(A2 — Ai)(l — /y)](i. Thus, r^j G [Airf, A2(i] where, we 
set the tube length d = 1 mm, Ai = 0.05, and A2 = 1.0. 
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FIG. 7. One of the generated IP patterns with g = 0.05 on 
a lattice of 200 x 300 nodes. The pattern was loaded into our 
network model. 



Above, rij is mapped to 1 — /^ because in the IP al- 
gorithm the next bond to be invaded is the one with the 
lowest threshold value, opposite to the network model, 
where the widest tubes will be invaded first. Note also, 
that in the network model the invasion of nonwetting 
fluid is controlled by the threshold capillary pressures pt 
of the tubes. According to Eq. (2) pt = 47/r in the 
middle of the tubes where x — d/2. In the IP model the 
distribution of /^ is flat. Thus, when r^ is mapped to /^ 
as described above, wc obtain a \/pt^ distribution of cap- 
illary pressure thresholds. However, since there is a one 
to one correspondence in the mapping between /^ and 
Pt, we can assume that the IP patterns are statistically 
equal to similar structures that would have been gener- 
ated in a full displacement simulation. The assumption 
provides that the displacement simulation is performed 
with an appropriate injection rate Q, according to g that 
was used to generate the IP patterns. 

After the IP patterns were successfully loaded into the 
network model, we started the simulations and ran the 
displacement a limited number of time steps while AP^ii 
was recorded. The number of time steps were chosen 
such that the front menisci got sufficient time to adjust 
according to the viscous pressure set up by the injection 
rate. For all four structures we chose M — 100 and Q — 
0.1 ml/min, giving Ca = 9.5x10"^. This Ca might be 
too high compared to the front widths we obtained at low 
Ca from simulations listed in Tables I and II. The reason 
why we choose a high Ca is to minimize computation 
time. Simulations show that fewer time steps and hence, 
less CPU time are required to adjust the front menisci 
when a high injection rate is applied instead of a low 
one. Moreover, the simulations also show that as long 
as the number of time steps are chosen sufficiently large 
to allow the front menisci to adjust, the exponent k in 
APj,|| oc A/i", is not sensitive on the injection rate. In the 
present simulations the number of time steps was 400. 

The result of the simulations is shown in Fig. 8 where 
we have plotted logio(APc||) versus \ogiQ{Ah). As for 
the previous results, wc find k = 1.0 ± 0.1. The slope 
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FIG. 8. logjQ(APcii) as function of logj^Q(A/i) for simula- 
tions initiated on IP patterns on lattices of 200 x 300 nodes. 
Ca = 9.5 X 10^'' and M = 100. The result is averaged over 
four different runs and the error bars denote the standard 
error in the mean. The slope of the straight solid line is 1.0. 



of the straight line in Fig. 8 is 1.0. We have also done 
displacement simulations on one of the IP patterns at 
Ca = 2 X 10~s with M = 1 and M = 100. These simula- 
tions were run in 1600 time steps and the result of those 
is consistent with Fig. 8. 



IV. EFFECT OF LOOPLESS STRANDS 

In [12] it was argued that AP^n — APnw — AP^ (see 
Fig. 1). At low injection rates or when the nonwet- 
ting phase is much more viscous than the wetting phase, 
APu. -C APnw, giving APc|| ~ APnw Thus, if the result 
of Xu et al. [12] should be valid for our problem, we would 
expect to find APj,|| oc Ah'^ where k = t/v + d^ — l — P/v. 
Inserting values of the exponents in 2D (t = 1.3, v = 4/3, 
d-E = 2, f3 = 5/36) gives k ~ 1.9. Our simulations clearly 
indicate that k ~ 1.0 which is inconsistent with the pro- 
posed result in [12]. Below we present an alternative view 
on the displacement pattern from that being initiated by 
Wilkinson [9] and used by Xu et al.. The alternative view 
is based upon the observation that nonwetting fluid flows 
in separate strands. 

Fig. 9 shows two typical displacement structures that 
were obtained from simulations at low and high Ca on the 
lattice of 40 x 60 nodes with M = 1 (Table III). We ob- 
serve that the nonwetting fluid (dark grey and black) gen- 
erates patterns containing no closed loops. That means, 
following a path on nonwetting fluid will never bring us 
back to the starting point. The loopless structure is a 
direct consequence of the evidence that a tube filled with 
wetting fluid and surrounded on both sides by nonwetting 
fluid is trapped due to volume conservation of wetting 
fluid. Because of trapped wetting fluid, the nonwetting 
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Ca = 1.6 X 10-5, M = 1.0 
FIG. 9. Two displacement structures of simulations at high 

Ca = 3.9x10"* (above) and low Ca = 1.6 xlO"'"' (below) be- 
fore breakthrough of nonwetting fluid. The nonwetting fluid 
(dark grey and black) is injected from below and wetting fluid 
(light grey) flows out along the top row. The lattice size was 
40 X 60 nodes and M = 1 (Table III). The black tubes de- 
note the loopless strands where nonwetting fluid flows and the 
dark grey tubes indicate nonwetting fluid unable to flow due 
to trapped regions of wetting fluid. Because of fluid trapping 
and subsequent volume conservation of wetting fluid, strands 
of different starting points along the inlet can never connect. 
Note the few fluid supplying strands from the inlet to the 
frontal region at low Ca compared to the case at high Ca- 



fluid also flows in separate strands, indicated as black 
tubes in Fig. 9. When the nonwetting fluid percolates 
the system there exists only on such strand connecting 
the inlet to the outlet. The dark grey tubes connecting 
to the strands are dead ends where nonwetting fluid can- 
not flow because of trapped wetting fluid. We note that 
the evidence of trapped wetting fluid in single tubes may 
easily be generalized to 3D and therefore our arguments 
should be valid there too. Similar loopless structures as 
in Fig. 9, were also pointed out in [27] for site-bond IP 
with trapping and in [28] for a loopless IP algorithm. 

From Fig. 9 we may separate the displacement pat- 
terns into two parts. One consisting of the frontal region 
continuously covering new tubes, and the other consist- 
ing of the more static structure behind the front. The 
frontal region is supplied by nonwetting fluid through a 
set of strands that connect the frontal region to the inlet. 
When the strands approach the frontal region they are 
more likely to split. Since we are dealing with a square 
lattice, a splitting strand may create either two or three 
new strands. As the strands proceed upwards in Fig. 9, 
repeatedly splits cause the frontal region to be completely 
covered by nonwetting strands. 

On IP patterns with trapping [27] or without 
loops [28,29] the length I of the minimum path between 
two points separated an Euclidean distance R scales like 
I oc R^' where Ds is the fractal dimension of the short- 
est path. We assume that the displacement pattern of 
the frontal region for length less than the correlation 
length (in our case Ws) is statistically equal to IP patterns 
in [27]. Therefore, the length of the nonwetting strands 
in the frontal region, is proportional to A/i^= where A/i 
is some vertical length less than Ws ■ If we assume that on 
the average every tube in the lattice has same mobility 
(kij/fiij), we obtain that the fluid pressure within one 
strand must drop like A/i" where k ~ Ds- Let us now 
consider the effect on the pressure when strands split. If 
we assume that the strands are straight {Dg = 1) then 
following a path where strands splits would cause the 
pressure to drop as Ah'^ where k < 1. This because the 
volume fluxes through the strands after a split must be 
less than the flux in the strand before it splits, due to 
volume conservation of nonwetting fluid. 

The two effects {k = Dg and k < 1) predict that 
the pressure drop in the nonwetting phase of the frontal 
region, AP„^, should scale as APnw oc A/i" where 
K < Dg- In 2D two different values for Dg have been 
reported: Dg = 1.22 [28,29] for loopless IP patterns, and 
Dg = 1.14 [27] for the single strand connecting the inlet 
to the outlet when nonwetting fluid percolates the sys- 
tem. We note that the result in [27] is essential equal 
to Dmin — 1-13 [25], that is the fractal dimension of the 
minimum path in 2D percolation where loops generally 
occur. Any of the above values for Dg together with 
the argument k < Dg, are supported by our simulations 
flnding k = 1.0 ± 0.1. 

Note the different pattern of strands at high and low 
Ca in Fig. 9. At low Ca few strands are supplying the 



frontal region with nonwetting fluid, and the strands split 
many times before the whole front is covered. At high 
Ca the horizontal distance between each strand in the 
static structure is much shorter, and only a few splits 
are required to cover the front. Moreover, we observe 
that at high Ca the length of individual strands in the 
front approaches the minimum length due to the tubes. 
In this limit we may treat the strands in the front as 
straight lines (i.e. Dg ^ \) causing k < 1. This is indeed 
supported by numerical results, finding that n decreases 
from about 1.0 to 0.8 when increasing Ca (see Fig. 2). 

Another important issue, arising at low Ca, is the effect 
of bursts on the capillary pressure. A burst occurs when 
a meniscus along the front becomes unstable and non- 
wetting fluid abruptly covers new tubes [22]. The strand 
where the burst initiates will during the burst, experi- 
ences a much higher fluid transport relative to strands 
far away. Describing the pressure behavior between the 
strand of the burst and the rest of the front is nontrivial. 
However, simulations show that even during bursts, we 
find that ^Pc\\ increases linearly with A/i. 

The indication that k ~ 1.0, may influence the scaling 
behavior of Ws as function of Ca- Assuming Darcy flow 
where the pressure drop depends linearly on the injection 
rate, we conjecture that APc|| fx CaAft.'*. Here APc|| de- 
notes the capillary pressure difference over a height Aft, 
when the front is stationary. That means, ^Pc\\ excludes 
situations where nonwetting fluid rapidly invades new 
tubes due to local instabilities (i.e. bursts). The above 
conjecture is supported by simulations showing that in 
the low Ca regime AP^i ^ Ca^h'^ where k ~ 1.0. Note, 

that APc|| 9^ ^-Pcjl in Fig. 2, since the latter includes 
both stable situations and bursts. 

At sufficiently low Ca the displacement may be mapped 
to percolation giving APc|| on f — fc oc S,^^^" [16,9,14]. 
Here / is the occupation probability of the bonds, fc is 
the percolation threshold, and ^ ex Wg is the correlation 
length. By combining the above relations for AP^n we 
obtain Wg oc Ca~°' where a = j^/(1 -t- vk). In 2D v — 4/3 
and inserting k = 1.0 gives a ~ 0.57. 

In Sec. Ill A we found that at high Ca the nonwetting 
fluid invades simultaneously everywhere along the front. 
Hence, the front never reaches a stationary state because 
of rapidly succeeding local instabilities. This is supported 
by simulations showing a crossover in APcii to a nonlinear 
dependency on Ca- Consequently, the above mapping 
to percolation might no longer be valid and we expect 
another type of functional behavior between Wg and C'a 
in the high Ca regime. 



V. COMPARISON WITH EXPERIMENTS 

Frette et al- [26] performed two phase drainage dis- 
placement experiments in a 2D porous medium with vis- 
cosity matched fluids (M = 1). They reported on the sta- 



• # Experiments 

• <> Simulations 




-5 -4 

log.o(C.) 

FIG. 10. logj^Q(ws) as function of logj,j(Ca) for experiments 
from [26] and simulations on the lattice of 40 x 60 nodes (Ta- 
ble III). For both experiments and simulations M = 1. The 
slope of the solid and dashed line is -0.6 and -0.3, respectively. 



bilization of the front and measured the saturated front 
width Wa , as function of Ca - For all our simulations ex- 
cept those performed on the IP patterns, wc have calcu- 
lated Wa- In Fig. 10 wc have plotted Wa as function of 
Ca in a logarithmic plot for the simulations in Table HI, 
(open diamonds) together with the experimental data of 
Frette et al (filled circles). 

In [26], their best estimate of the exponent a when 
assuming a power law Ws oc Ca^"^ was a ~ 0.6 ± 0.2, 
indicated by the solid line in Fig. 10. This is consistent 
with the suggested value a = i^/(l + vn) ~ 0.57 from 
Sec. IV. The simulations show a different behavior and 
they seem to fit a = 0.3 ± 0.1, according to the dashed 
line in Fig. 10. The simulations performed on the lattices 
of 25 X 35 nodes (Tables I and II) also give a ~ 0.3. 

Even though the overlap between experimental and nu- 
merical data in Fig. 10 is poor we suggest that the dif- 
ferent behavior of the experiments (at Ca ^ 1.0x10^^) 
and simulations (at Ca ^ 1.0x10"^) might be due to 
an expected change in a at high Ca- According to the 
discussion in Sec. IV it is not clear if the percolation 
approach giving a = v/{l + vn), is vahd for high Ca- 
The different scaling behavior observed in Fig. 10 might 
also be caused by the small system size of the simula- 
tions. At Ca — 1.0x10^"'' numerical simulations show 
that the front width becomes bounded by the system 
size, and therefore we are not able to observe a possible 
a ~ t^/(l + J^k) regime. We stress that more simula- 
tions on larger systems and at lower Ca are required in 
order obtain better overlap between simulations and ex- 
periments in Fig. 10. Until then, it is hard to draw any 
conclusions on the correct a. 

As a side mark, we note that our simulations giv- 
ing a ~ 0.3, are in agreement with numerical work 
in [12]. Their calculations oiwa were done for Ca between 
10^^ and 10^* coinciding with our region of simulations 
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in Fig. 10. According to Wilkinson [9] a = v/{l+t—l3+v) 
and by inserting values of the exponents in 2D we obtain 
a ~ 0.38. This is also within the uncertainties of our sim- 
ulation results. However, we emphasize that this might 
as well be a coincidence rather than an evidence, because 
Wilkinson's theory does not take into account that non- 
wetting fluid flows in strands along the front. 

A somewhat different process, but very interesting re- 
sult, is presented by Shaw in [30] . He measured the width 
of the drying front in a quasi 2D porous system and found 
that Ws oc i;^0-48±o.i jjgj-g yj jg lY^Q average front veloc- 
ity. Quite recently, this has been compared to theory 
in [31]. 



VI. CONCLUSION 

We have reported on the stabilization mechanisms of 
the front in drainage displacement going from low to high 
injection rates. The stabilization process was studied by 
using a network model simulating the viscous and cap- 
illary pressure buildup in the fluids during the displace- 
ments. We have found that the capillary pressure dif- 
ference AP(,|| 7 along the front varies almost linearly with 
the distance A/i, in the direction of the displacement. 



We conclude from simulations that AP, 



c|| 



/S.h'^ where 



our best estimate isK=1.0±0.1. This result supports 
the arguments showing n < Dg, where Dg is the fractal 
dimension of the loopless strands characterizing the dis- 
placement pattern. The evidence that nonwetting fluid 
flows in loopless strands along the front are not consid- 
ered in earlier proposed theories [9-12]. Hence, we con- 
clude that they are not compatible with drainage when 
nonwetting strands dominate the displacement process. 

Using the evidence that k ~ 1.0, we conjecture that 
the scaling of the front width Wg as function of Ca might 
alters from earlier suggestions in [9,11,12]. By mapping 
our problem to percolation we find Wg oc Ca""^ where 
a = V / {1 + vk). The result is consistent with experiments 
performed by Frette et al. [26] . Unfortunately, due to the 
small system sizes we are not able to confirm this scaling 
behavior by our simulations. We emphasize that a more 
stringent test on a should include simulations on larger 
systems and lower Ca, than presented here. 

In addition to AP^i we have calculated the capillary 
pressure variations along the front in the direction paral- 
lel to the inlet, APc±- Qualitatively, we have shown that 
APc± is a good indicator on whether the capillary pres- 
sures of the menisci along the front are all equal (capil- 
lary equilibrium) or fluctuating due to the viscous forces. 
When the capillary fluctuations are strong, we do not 
expect percolation to be a proper model for the displace- 
ment process. 
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APPENDIX A: 

Below we show how to deduce a — v/{l + t — j3 + v) 
in Wg oc Ca"" and find the corresponding exponent 
K — X+t/v+P/v in the power law APgii oc A/i'* when not 
considering that nonwetting fluid flows through strands. 
The calculations are carried out in two dimension, how- 
ever the extension to three dimensions is straight forward. 

Let us consider a piece of the nonwetting phase of size 
A/i in the frontal region. We assume that APc|| vary as 



AP^II oc wA/i'*, 



(Al) 



where v is the average fluid velocity in the pores. More- 
over, we assume that the front has reached a steady state 
and that the structure of the front is statistically equal 
to the front of an invasion percolation pattern. This as- 
sumption provides that A/i is sufficiently large for the 
percolation concept to apply but less than the front width 

Wg. 

The average nonwetting pore fluid velocity u, in the 
the region of size Aft,, is given by Darcy's law 



1 k APh 



S fi Ah 



(A2) 



Here S is the saturation of nonwetting phase, that is 
the volume fraction where nonwetting fluid can flow, and 
k is the permeability of the frontal region. According 
to percolation the frontal region is fractal, with fractal 
dimension D = d — P /v, giving 



and 



So,^^^—^ = Ah-^/-^, 



k oc A/i-*/". 



(A3) 



(A4) 



Here t is the conductivity exponent, /3 is the order param- 
eter exponent, and ly is the correlation length exponent 
in percolation. 

By inserting the expressions for S, fc, and AP^n into 
Eq. (A2) we find the exponent k = 1 + t/v — (3 /v. The 
exponent a follows by setting Ah = Wg and replace APcji 
in Eq. (Al) with the power law w,, oc ^ oc AP^'TJ''. Here ^ 
denote the correlation length in percolation. 
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